To APPEAR IN ApJ Letters 

Preprint typeset using 1^1^^ style cmulatcapj v. 4/9/03 



THE ORIGIN OF EPISODIC ACCRETION BURSTS IN THE EARLY STAGES OF STAR FORMATION 

E. I. VOROBYOV^'^, Shantanu Basu^ 
To appear in ApJ Letters 

ABSTRACT 

We study numerically the evolution of rotating cloud cores, from the collapse of a magnetically su- 
percritical core to the formation of a protostar and the development of a protostellar disk during the 
main accretion phase. We find that the disk quickly becomes unstable to the development of a spiral 
structure similar to that observed recently in AB Aurigae. A continuous infall of matter from the 
protostellar envelope makes the protostellar disk unstable, leading to spiral arms and the formation 
of dense protostellar/protoplanetary clumps within them. The growing strength of spiral arms and 
ensuing redistribution of mass and angular momentum creates a strong centrifugal disbalance in the 
disk and triggers bursts of mass accretion during which the dense protostellar/protoplanetary clumps 
fall onto the central protostar. These episodes of clump infall may manifest themselves as episodes 
of vigorous accretion rate (> 10^'' Mq yr^^) as is observed in FU Orionis variables. Between these 
accretion bursts, the protostar is characterized by a low accretion rate (< 10~^ Mq yr~^). During the 
phase of episodic accretion, the mass of the protostellar disk remains less than or comparable to the 
mass of the protostar. 

Subject headings: accretion, accretion disks — hydrodynamics — instabilities — ISM : clouds — MHD 
— stars: formation 



1. INTRODUCTION 

In our present view of low-mass star formation, a pro- 
tostar in the early stages of mass assembly (in the so- 
called class and class I phases) is surrounded by a 
[ protostellar disk which is in turn deeply embedded in 
an infalling envelope left over from the collapse of a 
rotating prestellar cloud core. The observed low lumi- 
nosity of these protostars implies a low mass accretion 
rate, and hence a long lifetime in order to achieve typ- 
ical final stellar masses; however, this is inconsistent 
with the number of known class and class I objects 
l)Kenvon et al.llT99(T[) . A possible explanation is that the 
protostellar accretion procee ds in two co-existing phases 
ijKenvon fc Hartmannlll995|) . Accretion from the enve- 
lope onto the protostellar disk takes place in a fairly 
uniform (though generally declining in time) manner, 
whereas accretion from the protostellar disk onto the 
central protostar occurs primarily in short (and infre- 
quently observed) but powerful episodes during which 
0.01 — 0.1 Mq can be accreted. These episodes of vigor- 
ous accretion M > 10""* Mq yr~^ manifest themselves as 
FU Orionis variables (FU Ori) . Between these accretion 
bursts, a typical class 0/class I protostar is characterized 
by a low accretion rate M ^ 10~^ Mq yr"-"^. 

The nature of FU Ori disk accretion bursts has been 
widely debated. For instance, close encounters in binary 
systems may cause a strong perturbation in protostellar 
disks and dri ye high accretion rates during a relatively 
short period f Bonnell fc Bast ion 1992). This mechanism 
requires rather eccentric orbits of binary systems and ob- 
viously fails to explain FU Ori outbursts in isolated pro- 
tostars. An alternative idea is that the thermal insta- 
bility (namely, the steep dependence of the disk opacity 
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on temperature between ~ 3000 K and ^ 10^ K) of the 
optically thick innermost regions of circumst ellar disks 
triggers FU Ori eruption s (Bell & Lin 1994: Cl arke et alJ 
ll99QtlLin fc Panaloizoulll985.) . Mos t thermal instability 
models exploit the a-prescription of lShakura &: SunvaevI 
()1973^ , who suggested that the disk effective viscosity is 
proportional to its temperature; an increase in disk tem- 
perature causes a higher rate of mass accretion due to 
an elevated viscous mass transfer and vice versa. Un- 
fortunately, many aspects of this mechanism of FU Ori 
outbursts are completely dependent upon the unknown 
value of disk effective viscosity, which determines the 
timescale for outbursts. 

It is known from theoretical and numerical stud- 
ies that protostellar disks may be subject to the de- 
velopment of global spiral instabilities. For instance, 
Eaughlin & Bodcnhcimcr ( 1994) have studied the nonax- 
isymmetric evolution of protostellar disks and found that 
they are susceptible to a series of spiral instabilities, with 
the fastest growing modes being the one-armed (m = 1) 
and two-armed (m = 2) patterns. Recent numerica l 
simulations of star cluster formation l)Bate et al.l 120031) 
confirmed that the protostellar disks formed around pro- 
tostars were gravitationally unstable and prone to the 
development of spiral density waves. Simulations of 
marginally unstable protoplanetary disks show the for- 
mation of flocculent and clumpy spiral structure, sug- 
gesting a possible transient rise in the mass a ccretion rate 
associated with clump infall (,Boss .2003l) . iMeiia et alJ 
(2005]) have also reported three-dimensional simulations 
which demonstrate a single FU-Ori-like outburst asso- 
ciated with the growth of spiral structure in an iso- 
lated protoplanetary disk. Recent observations do reveal 
a flocculent spiral structu re in the proto stellar disk of 
AB Aurigae l|Corder et al.l l2004: Fukaga wa et alJ l2004'l. 

In this Letter, we present the first model of cloud core 
collapse which self-consistently generates multiple accre- 
tion bursts that can be identified with FU Ori eruptions. 
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The protostellar disk in our model is formed as a result of 
the collapse and is not isolated from the parent core en- 
velope. The details of the numerical model are described 
in § 121 The results of simulations are presented in § 
Our main conclusions are summarized in § 01 

2. MODEL DESCRIPTION 

We model the collapse of a rotating cloud core which 
is threaded by a frozen-in magnetic field with spatially 
uniform mass-to-flux ratio. The magnetic field effect is 
comparable to but weaker than gravity, so that the core is 
magnetically supercritical. A magnetically-diluted grav- 
itational collapse ensues. Our initial core model is a 
good approximation to the supercritical cores that result 
from the fragmentation of magnetized clouds that are 
initially either critical/s ubcritical or supercritical (e.g., 
[Basu f997; Ba su fc Ciol ek 2004). 

We follow the collapse through the prestellar collapse 
phase and into the accretion phase which sees the de- 
velopment of a protostar and a protostellar disk. The 
thin-disk approximation is used in the form appropriate 
for a supercritical co re with ma s s-to-fiux ratio m that 
is spatially uniform jBasulll997t iNakamura & Hanawal 
119971 ISbn fc"illT99l . The magnetic field pressure en- 
hances the gas pressure and the magnetic tension dilutes 
the effect of gravity. The disk is symmetric about the 
midplane, with a vertical magnetic field component 
inside the disk and both vertical and tangential compo- 
nents outside it. The material external to the disk is 
current-free (j — 0). We characterize the field strength 
by the parameter a = fi^^ — -Bz/(27rVGS), where S is 
the gas surface density; note a < f for a supercritical 
core. 

The magnetohydrodynamic equations we use are writ- 
ten in the thin-disk approximation as 
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where v is the gas velocity in the disk plane, P is the 
vertically integrated gas pressure, and $ is the gravita- 
tional potential. Equation ^ contains the Lagrangian 
derivative d/dt — d/dt + v ■ W. 

We assume a two-component equation of state P = 
CsS-|-CsScr(S/Scr)''', whcrc Cs is the isothermal speed of 
sound and Scr is the critical gas surface density above 
which the disk becomes optically thick. This equation of 
state allows for a smooth transition between the isother- 
mal and non-isothermal regimes during the collapse. We 
use a canonical val ue o f the critical gas volume density 
ricr — 10^^ cm~^ ((Larsonll2003il . which is equivalent to 
Ecr = 36.2 g cm~^ for the gas disk in vertical hydrostatic 
equilibrium. In equation (0), we assume a cloud scale 
height Z = c2/(7rG'S) for S < Scr and Z = c^/iirGl^cr) 
for E > Ecr- We adopt the ratio of specific heats 7 = 7/5 
for the optically thick regime, appropriate for an adia- 
batic diatomic gas. 

Equations and (O are numerically solved in po- 
lar coordinates (r, 0) using the method of finite differ- 
ences with a time-explicit, operator-split solution proce- 
dure similar to that d escribed by Sto ne and Norman in 
their ZEUS-2D code l|Stone fc Norman 1992) . The de- 
tails of the code and results of the tests will be given in a 



follow-up paper. The numerical grid has 256 x 256 points, 
which are logarithmically spaced in r-direction allowing 
for a good resolution in the inner cloud regions. The in- 
nermost grid point is located at 10 AU and the size of 
the first adjace nt cell is 0.3 AU. The Truelove criterion 
l)Trnelove et al.t.1998'1 is preserved throughout the simu- 
lations - the size of grid cells is smaller than the Jeans 
length. We introduce a "sink cell" at r < 10 AU, which 
represents the central protostar plus some circumstellar 
disk material, and impose a free inflow inner boundary 
condition. We assume that the matter is cycled through 
the circumstellar disk and onto the protostar rapidly 
enough so that the mass infall through the sink cell is 
at least proportional to the mass accretion rate onto the 
protostar. We impose the outer boundary condition such 
that the gravitationally bound cloud has a constant mass 
and volume. The gravitational potential of the cloud is 
evaluated numerically usin g the fast Fourier transform 
l|Binnev fc Tremainelll987|l . 

3. RESULTS 

We have studied many different initial cloud configura- 
tions and present here a prototype magnetized (a = 0.45) 
rotating cloud with mass M^i = 2.45 Mq and an outer ra- 
dius Tout — 20000 AU. The cloud is composed of molecu- 
lar hydrogen with a 10% admixture of atomic helium and 
is initially isothermal at T = 10 K (cs = 0.188 km s^^). 
The initial surface density and angular velocity distribu- 
tions are those chara cteristic of a collapsing axisymmet- 
ric supercritical core llBasulll997D : 
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Here, Eq = 3.55 x 10~ g cm^"^ and fio = 0.5 km s~^ pc~ 
are the central surface density and angular velocity, re- 
spectively. We choose a value tq = Cg/{1.5GT,o), so that 
ro is comparable to the Jeans length of an isothermal 
sheet. We mimic the slight non axisymmetry in more re- 
alistic models of core formation ("Bas u & Ciolekll2004H by 
substituting in Eq. © with r^(cos^ -|- sin^ (j)), 
where the parameter a = 0.98 denotes the cloud oblate- 
ness. The ratios of rotational, magnetic, and thermal en- 
ergies to the gravitational energy of the cloud are 0.4%, 
30%, and 27.2%, respectively. Thus, the initial cloud is 
gravitationally unstable. We emphasize that our quali- 
tative results are insensitive to the particular choice of 
initial conditions. 

The cloud evolution is characterized by a slow ini- 
tial gravitational contraction and then a very rapid run- 
away collapse until the formation of the central protostar. 
The black line in Figure ^ shows two distinct phases 
in the temporal evolution of the mass accretion rate M 
onto the protostar. The early behavior of M is qual- 
itatively similar to t hat obtained in spherica l collapse 
simulations (see e.g., iVorobvov fc BasiJ I2005D . Accre- 
tion shows a very rapid increase to a maximum value 
of M = 1.0 X 10-* Mq yr-i at t = yr, when the 
central protostar forms. Subsequently, there is a slow 
decline in M, when the gas is accreted directly onto the 
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Fig. 1. — Periodic mass accretion outbursts. The temporal evolution of a) the mass accretion rate M (the black line) and the Toomre 
parameter Q (the red line) and b) the bolometric luminosity (the brown line) and the normalized gravitational torque per unit mass 
r/r max (the blue line). The horizontal axis shows the elapsed time since the formation of the protostar. The behavior of M shows two 
distinct phases. In the early phase, M slowly declines and tends to approach a constant value. In the later phase, after the formation of 
the protostellar disk at t = 7700 yr, the mass accretion occurs in periodic bursts. The evolution of both Q and r/Fmax show a correlation 
with the accretion bursts. 



protostar from the inner envelope, which has a relatively 
low specific angular momentum. The second phase starts 
at t = 7700 yr when the protostellar disk forms around 
the protostar due to the infall of matter from the outer 
envelope, which has a higher specific angular momen- 
tum. The mass of the protostar at this stage is approx- 
imately Ms — 0.4 Mq. The matter in the disk moves 
on nearly circular orbits and most of it is far from the 
protostar. Hence, the accretion rate onto the protostar 
abruptly drops down to a negligible value. A continuous 
infall of matter from the envelope makes the protostel- 
lar disk unstable to the development of spiral structure 
shown in Fig. [21 and induces the formation of dense pro- 
tostellar/protoplanetary clumps within the arms. The 
even (m = 2,4) modes are dominant in this simula- 



tion, but the odd (m ==1,3) modes are also excited in 
some simulations. Spiral arms tr ansport an gular m omen- 
tum outward and mass inward ifLvnden-Bell fc KalnaisI 
Il972j) . The ensuing redistribution of mass and angu- 
lar momentum creates a strong centrifugal disbalance 
in the protostellar disk and triggers bursts of mass ac- 
cretion when dense protostellar/protoplanetary clumps 
in the inner disk are driven into the protostar (an 
animation of the disk evolution can be downloaded 
from http : //www . astro . uwo . ca/^basu/mv . htm). Dur- 
ing this process, the mass in the protostellar disk re- 
mains somewhat less than the mass of the protostar. 
The episodes of clump infall can manifest themselves as 
very short (< 100 yr) but vigorous {M = [1 — 10] x 
10~^ Mq yr~^) accretion bursts as is clearly seen in 
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Fig. 2. — A sequence of gas surface density images showing the evolution of the protostellar disk after the formation of the protostar at 
t = yr. The disk quickly becomes unstable and develops a spiral structure with dominant m = 2 and m = 4 modes. Formation of dense 
protostellar /protoplanetary clumps within the spiral arms is evident in most images. The numbers in the left upper corner of each image 
show the elapsed time since the formation of the protostar. 



Fig. n During the accretion bursts, 0.01 — 0.05 Mq of 
gas is accreted and the accretion himinosity may by grow 
many orders of magnitude. The duration of the inter- 
vening quiescent accretion phase with M = (1 — 10) x 
10~^ Mq yr~^ is usually (1 — 3) x 10"^ yr. The frequency of 
bursts decreases with time and the number of bursts may 
amount to 15-30. The brown line in Fig.^ shows a lumi- 
nosity Lboi = GMgM /Rc (it is assumed to derive entirely 
from the disk ac cretion), where Rr = ^R(^ is the radius 
of the protostar l)Masunaga fc lnutsukall20do|) . Our Lboi 
is an upper limit to the expected observable bolometric 
luminosity. It is evident that the episodes of clump in- 



fall lead to a dramatic increase in Lboi (by a factor of 
up to 2000 as compared to a quiescent period) and the 
protostar may reveal itself as an FU Orionis variable. 

Our calculations o f the approxim ate Toomre parame- 
ter Q = Csf^/(7rGE) (jToomrdll98l|) and the total gravi- 
tational torque per unit mass F, which is the sum of the 
individual torques per unit mass {\d^/d4>\) on all com- 
putational cells, support the above scenario of episodic 
accretion. The quantity Cg = dP/dT, is the squared ef- 
fective sound speed. The Q parameter may serve as an 
approximate stability criterion - gas disks are gravita- 
tionally unstable to local nonaxisymmetric perturbations 
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if Q < 1.5-1.7 llBosslll998tlNelson et aLlllflOSl) . while T 
may roughly express the efficiency of angular momen- 
tum and mass redistribution by spiral inhomogeneities 
in the disk (Tomlcv ct al. 1994). The Toomre parameter 
is calculated by averaging Cs, f^, and S over all compu- 
tational cells. The red and blue lines in Fig. QJl and^^ 
show the evolution of Q and the normalized gravitational 
torque F/Fmax after the protostellar disk formation, re- 
spectively. In the early phase of near constant accretion, 
the matter is directly accreted onto the protostar and Q 
is much larger than unity. When the protostellar disk 
forms, its density starts to grow due to accretion. As a 
consequence, the Toomre parameter gradually decreases 
below the stability limit Q ~ 1.5 and reaches a minimum 
value at the time of the accretion burst. This strongly 
suggests a causal link between the gravitational instabil- 
ities and accretion bursts. The behavior of F also shows 
a direct correlation with accretion bursts. The gravita- 
tional torque gradually increases and reaches a maximum 
value at the time of each accretion burst, indicating the 
growing efficiency of inward mass transport before the 
burst. We note that the strength of the torque due to 
artificial viscosity (used in our Eulerian code to smooth 
shocks) is at least an order of magnitude smaller than the 
strength of the gravitational torque associated with spi- 
ral instabilities. Thus, the artificial viscosity cannot be 
responsible for the accretion bursts. After the burst, the 
mass of the protostellar disk decreases and the growth 
of spiral instabilities becomes temporarily suppressed, as 
indicated by high values of Q (> 1.5) and a sharp de- 
crease in F. However, the continuous mass infall onto 
the disk from the envelope quickly destabilizes the disk 
and the cycle repeats until most of the envelope mass is 
accreted by the protostar. 

Ambipolar diffusion, while not included in our model, 
is expected to favor the formation of clumps and sub- 
sequent burst activity by removing magnetic support. 
Indeed, a model with no magnetic support {a = 0) 
shows an increase in the frequency and amplitude 
of accretion bursts (Vorobyov &: Basu, in prepara- 
tion) . Further magnetic effects such as magnetic braking 
ijKrasn opolskv fc Konigl' '2002^ and magnetorotational 
instability (,Fromang et al 2005) may become important 



in the late accretion phase, but can only be studied in a 
future three-dimensional model. 

4. CONCLUSIONS 

Rotating protostellar cores show two distinct phases 
in the temporal evolution of the mass accretion rate M 
onto the protostar. The early behavior of M is qual- 
itatively similar to that obtained in spherical collapse 
simulations. Accretion shows a very rapid increase to a 
maximum, when the central protostar forms, and a sub- 
sequent slow decline, when the gas is accreted directly 
onto the protostar from the inner envelope. The second 
phase starts when the protostellar disk forms around the 
protostar due to the infall of matter from the outer en- 
velope with higher specific angular momentum. In this 
phase, M is characterized by very short (< 100 yr) but 
vigorous (M = [1 — 10] x 10"'' Mq yr~^) accretion bursts, 
which are intervened with longer periods (~ 10'^ yr) of 
quiescent accretion. 

We have demonstrated that the repeating accretion 
bursts reflect a basic self- regulation mechanism that is in- 
herent to self-gravitating rotating protostellar disks. We 
emphasize that it is the ongoing infall of matter from 
the protostellar envelope that continually destabilizes 
the disk and causes it to periodically dump significant 
amounts of matter onto the protostar while transferring 
excess angular momentum to the envelope. The effect of 
additional support due to a frozen-in supercritical mag- 
netic field is to moderate the burst activity but not sup- 
press it. Recent high resolution Hubble Space Telescope 
observations of the central regions (< a few h undred pc) 
of Seyfert galaxies ijR.egan fc Mulchaevlll999D also reveal 
rich spiral structures. We suggest that a self-regulation 
mechanism similar to that in our model may operate on 
galactic scales and be responsible for periodic nuclear ac- 
tivity in at least some Seyfert galaxies. 



This research was supported by the Natural Sciences 
and Engineering Research Council of Canada. EIV grate- 
fully acknowledges support from a CITA National Fel- 
lowship. 



REFERENCES 



Basu, S. 1997, ApJ, 485, 240 

Basu, S., & Ciolek, G. E. 2004, ApJ, 607, L39 

Bate, M. R., Bonnell, I. A., & Bromm, V. 2003, MNRAS, 339, 577 

Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987 

Binney, J., & Tremaine, S. 1987, Galactic Dynamics, (Princeton 

Univ. Press, Princeton), 96 
Bonnell, I., & Bastien, P. 1992, ApJ, 401, L31 
Boss, A. P. 1998, ApJ, 503, 923 
Boss, A. P. 2003, ApJ, 599, 577 

Clarke, C. J., Lin, D. N. C, & Pringle, J. E. 1990, MNRAS, 242, 
439 

Corder, S., Eisner, J., & Sargent, A. 2005, ApJ, 622, 133 
Fromang, S., Balbus, S. A., Terquem, C, & De Villiers, J.-P. 2005, 
ApJ, 616, 364 

Fukagawa, M., Hayashi, M., Tamura, M. et al. 2004, ApJ, 605, L53 
Kenyon, S. J., Hartmann, L. W., Strom, K. M., & Strom, S. E. 

1990, AJ, 99, 869 
Kenyon, S. J., & Hartmann, L. W. 1995, ApJS, 101, 117 
Krasnopolsky, R., & Konigl, A. 2002, ApJ, 580, 987 
Larson, R. B. 2003, Rep. Prog. Phys., 66, 1651 
Laughlin, G., & Bodenheimer, P. 1994, ApJ, 436, 335 



Lin, D. N. C. & Papaloizou, J. C. B. 1985, 
in Protostars and Planets II, ed. D. C. Black, M. C. Matthews, 
Eds. (Univ. Ariz. Press, Tucson, 1985), 981 

Lynden-Bell, D., & Kalnajs, A. J. 1972, MNRAS, 157, 1 

Masunaga, H., & Inutsuka, S. 2000, ApJ, 531, 350 

Mejia, A. C, Durisen, R. H., Pickett, M. K., & Cai, K. 2005, ApJ, 
619, 1098 

Nakamura, F., & Hanawa, T. 1997, ApJ, 480, 701 
Nelson, A. F., Benz, W., Adams, F. C, & Arnett, D. 1998, ApJ, 
502, 342 

Regan, M. W., & Mulchaey J. S. 1999, AJ, 117, 2676 
Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 
Shu, F. H., & Li, Z.-Y. 1997, ApJ, 475, 251 
Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753 
Tomley, L., Steiman-Cameron, T. Y., & Cassen, P. 1994, ApJ, 422, 
850 

Toomre, A. 1981, in The Structure and Evolution of Normal Galaxies, 
ed. S. M. Fall, D. Lynden-Bell, (Cambridge Univ. Press, 
Cambridge), 111 

Truelove, J. K., Klein, R. I., McKee, C. F., HoUiman II, J. H., 
Howell, L. H., Greenough, J. A., & Woods, D. T. 1998, ApJ, 
495, 821 

Vorobyov, E. I., & Basu, S. 2005, MNRAS, 360, 675 



